function yyy=f0_August25_joint_eight_parameters_calibration_function(xxx,params)

    alpha=params.alpha;
    pi=params.pi;
    eta=params.eta;
    r=params.r;
    
    params.lambdaE=xxx(1);
    params.lambdaS=xxx(2);
    params.lambdaI=xxx(3);
    params.lambdaC=xxx(4);
    params.f=0;
    params.xiE=xxx(5);
    params.xiI=xxx(6);
    params.lambdae=params.lambdaE;
    params.xie=xxx(7);
    
    guess=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx,fval] = fminsearchbnd(@(x) equation_difference(x,params),guess,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    alpha=params.alpha;
    pi=params.pi;
    eta=params.eta;
    r=params.r;
    lambdaE=params.lambdaE;
    lambdaS=params.lambdaS;
    lambdaI=params.lambdaI;
    lambdaC=params.lambdaC;
    lambdae=params.lambdae;
    f=params.f;
    xiE=params.xiE;
    xiI=params.xiI;
    xie=params.xie;
    
    xI=xx(2);
    xE=xx(5);
    xe=xx(8);
    tau=xE+xe;
    theta=xe/(xe+xE+xI);
    B=xx(6);
    gq=xx(7);
    RD_all=xiI*eta^(-1)*(xI)^(1/(1-alpha))+xiE*eta^(-1)*(xE)^(1/(1-alpha))+xie*eta^(-1)*(xe)^(1/(1-alpha))+B*xI*eta^(-1)*pi^(-1)*(1-f);
    RD_new=xie*eta^(-1)*(xe)^(1/(1-alpha));

    moment1=xI*(eta-1)*log(1-lambdaC)+(xE+xe)*(eta-1)*log(1-lambdaS); % target -0.1506
    moment2=((xE+xe)*log(1-lambdaS))/(xI*log(1-lambdaC)); % target 1.4858
    moment3=xE*((1+lambdaE)^(eta-1)+(1-lambdaS)^(eta-1)-2)+xI*((1+lambdaI)^(eta-1)+(1-lambdaC)^(eta-1)-2)+xe*((1+lambdae)^(eta-1)+(1-lambdaS)^(eta-1)-2); % target 0.016
    moment4=xe*(1+lambdae)^(eta-1); % target 0.0114
    moment5=xI/xE; % target 6
    moment6=RD_all; % target 0.0988
    moment7=xI+xE; % target 0.145
    moment8=xe; % target 0.0108
    moment9=RD_new; % target 0.0115

    weight=params.weight;
    %yyy=weight(1)*(abs(moment1+0.1506)/0.1506) + weight(2)*(abs(moment2-1.4858)/1.4858) + weight(3)*(abs(moment3-0.016)/0.016) + weight(4)*(abs(moment4-0.0114)/0.0114) + weight(5)*(abs(moment5-6)/6) + weight(6)*(abs(moment6-0.0988)/0.0988) + weight(7)*(abs(moment7-0.145)/0.145) + weight(8)*(abs(moment8-0.0108)/0.0108) + weight(9)*(abs(moment9-0.0115)/0.0115);
    yyy=weight(1)*(abs(moment1+0.1506)/0.1506) + weight(2)*(abs(moment2-1.4858)/1.4858) + weight(3)*(abs(moment3-0.016)/0.016) + weight(4)*(abs(moment4-0.0114)/0.0114) + weight(6)*(abs(moment6-0.0988)/0.0988) + weight(7)*(abs(moment7-0.145)/0.145) + weight(8)*(abs(moment8-0.0108)/0.0108) + weight(9)*(abs(moment9-0.0115)/0.0115);

    if fval>1e-15
        yyy=10^10; % panelize the case where equlibrium solver doesn't find equlibrium under certain paramters
    end
    